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We report Brownian dynamics (BD) simulation and theoretical results for a system of spherical colloidal particles with permanent 
dipole moments in a rotating magnetic field. Performing simulations at a fixed packing fraction and dipole coupling parameter, 
we construct a full non-equilibrium phase diagram as function of the driving frequency (wq) and field strength {Bq). This diagram 
contains both synchronized states, where the individual particles follow the field with (on average) constant phase difference, and 
asynchronous states. The synchronization is accompanied by layer formation, i.e. by spatial symmetry-breaking, similar to 
systems of induced dipoles in rotating fields. In the permanent-dipole case, however, too large loq yield a breakdown of layering, 
supplemented by complex changes of the single-particle rotational dynamics from synchronous to asynchronous behavior We 
show that the limit frequencies can be well described as a bifurcation in the nonlinear equation of motion of a single particle 
rotating in a viscous medium. Finally, we present a simple density functional theory, which describes the emergence of layers in 
perfectly synchronized states as an equilibrium phase transition. 



1 Introduction 

The dynamics of anisotropic particles driven by time- 
dependent, magnetic or electric, external fields is currently a 
I topic receiving much attention. Many experimental and theo- 
retical studies in this area focus on the field-induced dynamics 
■of a isolated nanoparticle such as a magnetic rod," a mag- 
netic chain"* or filament, ^ or an optically excitable nanorod- in 
a viscous medium. Understanding the resulting single-particle 
rotational dynamics is particularly important for actuators,^ 
molecular switches, particles in optical traps, and in the more 
general context of microfluidics. ^ From the theoretical side, 
these problems are often successfully analyzed on the basis of 
single-particle, nonlinear equations for the driven rotational 
motion in the presence of solvent-induced friction.— ii£ Typi- 
cally, the particle dynamics exhibits a "linear" regime at low 
driving frequencies, where the particle axis follows the field, 
and various types of nonlinear behavior at high frequencies, 
such as rotation against the torque. Many of these nonlinear 
phenomena, including transient behavior such as conformal 
transitions'* of magnetic chains following a sudden switch-on 
of the driving field, can also be observed experimentally.-i*^ 

Apart from the single-particle dynamics, another current 
focus concerns the self-assembly behavior in colloidal many- 
particle systems that are exposed to rotating fields. Indeed, in 
material science, time-dependent fields are currently realized 
as a powerful tool to control self-assembly processes, which 
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are an important prerequisite for synthesizing functional mate- 
rials.^-^ A classical example in this context, first discussed by 
Martin et al.,^ are systems of paramagnetic (or polarizable) 
spherical particles in magnetic (electric) fields rotating in a 
plane. For sufficient field strength, both experiment and com- 
puter simulations ^ ' ' reveal the formation of layers in the field 
plane, i.e. a spatial symmetry breaking induced by the rotating 
field. Indeed, a rotating in-plane field generates, on averaging 
over time, an inverted dipolar pair interaction with in-plane at- 
traction and repulsion along the rotation axis.Jiii^ Therefore, 
the structures induced by planar rotating fields are markedly 
different from those observed in a consant and homogeneous 
field, which supports the formation of field-aligned chains 
(low densities) '^'"'^ and bulk crystals.i^iii 

The general idea to use time-dependent fields to tune 
pair interactions and thereby control the morphology of self- 
assembled structures has meanwhile become more and more 
popular (see Refs. 18, 19), a recent example being the forma- 
tion of self-healing membranes of superparamagnetic particles 
in tilted rotating fields. Interestingly, these self-assembly 
phenomena can often be explained from an equilibrium per- 
spective involving the free energy and resulting phase behav- 
ior of a many-particle system in a time-averaged field. Clearly, 
the crucial assumption in adopting this perspective, which is 
also often exploited in computer simulation studies (see e.g. 
Refs. mil) is that the particles follow the field synchronously. 

In the present paper we explore, for a magnetic many- 
particle system, the link between the collective, self-assembly 
behavior, on the one hand, and the single-particle dynam- 
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ics, on the other hand. Specifically, we consider a ferrofluid 
subjected to a rotating in-plane field, where the ferrofluid 
is modeled by a system of dipolar soft spheres (DSS). The 
same model has been considered earlier in a computer sim- 
ulation study by Murashov and Patey,^^ where the aim was 
to demonstrate that layering occurs not only in systems of 
(super-)paramagnetic or polarizable particles (as those consid- 
ered by Martin et a/. '^''''), but also for particles with perma- 
nent dipoles. Here we investigate the driven DSS system both 
by Brownian dynamics (BD) computer simulations, which are 
described in Sec. IH and by theory. As a first main result, we 
present in Sec. 13. 21 a full non-equilibrium "phase" diagram in- 
dicating the domain of layer formation in the plane spanned 
by the frequency and strength of the driving field at constant 
equilibrium thermodynamic parameters. Secondly, to iden- 
tify the role of mutual synchronization of the particles, we in- 
vestigate in Sec. 13.31 the rotational dynamics within layered 
and unlayered states by analyzing suitable distribution func- 
tions. A similar strategy has recently been proposed in a dy- 
namic density functional study of rod-like particles in rotating 
fields.— In Sec. l3.3.ll we show that the breakdown of layering 
observed at high frequencies in the ferrofluid system can be 
described by a single-particle theory similar to those used for 
field-driven single nanoparticles in viscous media.— Finally, in 
Sec. 13.41 we propose a simple equilibrium density functional 
approach to investigate the role of translational entropy for 
layering in synchronized ferrofluid systems. The results are 
in good agreement with corresponding BD simulations. We 
close the paper with a brief summary and conclusions (Sec. 5). 

2 Model and simulation methods 

In our simulations we model the colloidal suspension by a sys- 
tem of dipolar soft spheres (DSS). The solvent is not explic- 
itly taken into account. The DSS pair potential between two 
spheres is comprised of a repulsive potential [/'°p and a point 
dipole-dipole interaction potential 

u-^in,)~'^'--^f--^^Kf^. (1) 

In eqn (HJ, r^j is the vector between the positions of the parti- 
cles i and j, r^j its absolute value, and fi^ is the dipole moment 
of the ith particle. The potential W^^p is the shifted soft sphere 
potential, which is given by 



C/-P(r) = U^^ir) - U^^ir,) + (r, ~ r)::^(r,), (2) 
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is the unshifted soft sphere (SS) potential for particles of ra- 
dius cr. Further, Vc — 2.5(7 is the radius at which we cut off 
the potential C/"=p. 

We investigate the system using non-overdamped Brownian 
dynamics (BD) simulations (sometimes called Langevin dy- 
namics simulations). The corresponding equations of motion 
for particles of mass m and moment of inertia / are^ 



Iu>i 



pDSS 

rpDSSJ 



(4) 
(5) 



In these equations — ksT/D and = ksT/Dr are 
friction coefficients with fc^ and T being Boltzmann's con- 
stant and temperature, respectively, while D and Dr are the 
translational and rotational diffusion constants. Furthermore 
uji is the angular velocity of particle i, Yf and Tf are ran- 
dom Gaussian forces and torques, T^'^* are torques due to an 
external field, and = — (/Xj • Ti)/Xj//^|. Their cartesian 
components (a, (3 — x,y, z) satisfy 



{Fg{t))=Q 
{TF^m = 



as well as 



{Fg{t)Ffp{t')) - 6fcBrCT<5„<5„/3<5(t - t') 
{Tg{t)Tf^At')) - QkBTiR5,,5^p6{t - t'). 



(6) 
(7) 



(8) 
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As eqns ©-([Oil show, the friction coefficients and the prob- 
ability distributions of the random forces and torques are re- 
lated via the fluctuation-dissipation theorem. This ensures that 
the system approaches a canonical distribution of states char- 
acterized by a constant temperature T in the absence of an 
external drive. To deal with the long-ranged dipolar interac- 
tions, we used the Ewald summation method with conduct- 
ing boundaries.^"* We have parallelized the evaluation of the 
Ewald sum in our simulation with OpenMP and MPI. The 
equations of motion were integrated with a Leapfrog algo- 
rithm. 

The external field that the particles interact with rotates with 
frequency ojo in the x — y-plane and is given by 



B [t) = Bq {e^ cos loqI + ey sin ujQt) . 



(10) 



For convenience, we make use of the following reduced 
units: Field strength Bq — (a^/eY^'^Bo; angular frequency 
ujq = {ma^ /ey^^LUo; density p* — p; dipole moment 
/i* = (ecr"^)^^/^^; and moment of inertia /* = {ma^)~^I. 
Unless stated otherwise, the simulations were carried out with 
864 particles at density p* — 0.1, dipole moment p* — 3, mo- 
ment of inertia /* = 0.025, and temperature T* = ksT/e = 
1.35. To verify our results, we also ran simulations with 
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up to 4000 particles. The translational and rotational diffu- 
sion constant were chosen to be — 0.1 • (ea^/m)^/^ and 
Dr = 3 ■ (m(T^/e)^^/^, respectively, and we used a timestep 
of At = 0.0025 • (mcr^/e)^/^. These values are consistent 
with those chosen in earlier BD studies of rotating dipolar 
systems.— We note, however, that the effects reported in the 
present paper also appear for other values of D and D^- 



3 Results and discussion 
3.1 Zero field system 

The zero field system, which represents our starting point, 
is characterized by a large dipolar coupling strength A = 
M^/ ikBTa^) ~ 6.7 and a relatively low density. As expected 
for such a strongly coupled system, the particles self-assemble 
into chainlike structures. '^-^^ This can be seen in the snap- 
shot depicted in Fig. [T^. Our reason for considering a sys- 
tem of a coupling strength this high is that this seems to be 
a prerequisite for layer formation. Indeed, irrespective of the 
field strength, we did not observe any layering for values of 
A that are smaller than approximately 4.6 (at the temperature 
T* = 1.35). 

Contrary to A, our choice of the density is less restricted, 
since the layering phenomenon persists over a wide range of 
densities (at least up to p* ~ 0.4). However, choosing the 
small density of p* — 0.1 has the advantage that layers are 
easily discernible. 




Fig. 1 (a) Snapshot of the system in zero field at = 0.1, 

T* — 1.35, and p* = 3. (b) Snapshot of a system in a layered state. 

The strength and frequency of the field are Bq ~ 12 and cja = 15, 

respectively. 



3.2 The layering effect 

We now consider the same system in rotating fields of various 
strengths B'q and frequencies Wq. For sufficiently large Bq 
and not too high frequencies (see below), the particles arrange 
themselves into layers. An example of this is shown in Fig. [TJi. 



This phenomenon was first explained by Halsey, Anderson, 
and Martin. They realized that the time-averaged potential 
between two particles i and j that rotate with the same angular 
frequency (given by the external field) and are aligned with 
each other, i.e. rotate circularly in a synchronized fashion with 

Hi{t) = Hj{t) — p{ex coswo^ + Gj/ sincjot), (H) 
is given by 



/ C/°(r,,,/x,(t),/x^-(i))di 

Jtn 

2(l-3cos2e,j) 
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(12) 



In this equation, is the dipole-dipole potential (see eqn 
([T]), T = 27r/ajo is the oscillation period, and is the angle 
between the interparticle vector Yij and the direction perpen- 
dicular to the plane of the field. As shown by the last line 
in eqn (fT2l l. the time-averaged potential corresponds to an in- 
verted dipolar (ID) potential, which is attractive if the angle 
Qij satisfies cos^ Qij < 1/3, i.e. if the particles i and j are 
approximately in the same plane with respect to the field. Con- 
versely, if the angle Qij satisfies 1/3 < cos^ 8^^, the particles 
repel each other This direction dependence of the ID potential 
explains why layers are a favorable configuration for a driven 
system in which the particles rotate synchronously. 

Note that for the above argument to hold, the translational 
motion of the particles should be small compared to their ro- 
tational motion. More precisely, during one rotational period 
they should migrate much less than their own diameter 

In the following we aim to determine more precisely the 
range of frequencies and field strengths at which layering oc- 
curs. To do that, we need a suitably defined order parameter 
We tested several ones and compared them with one another 
The order parameter that we will use here is given by 



1 ^ 



(13) 



i=l 



where is the total number of particles, ( • • • ) denotes a time- 
average, and rii is defined as follows: Consider a sphere of 
radius tq around particle i. Divide that sphere into two parts, 
one of which is given by the points within the sphere whose 
distance vector to particle i together with the z-axis encloses 
an angle Q satisfying —0.5 < cos0 < 0.5 (see Fig. |2|i. If 
there are more (less) particles in this equatorial volume than 
in the polar volume around particle i, set = 1 (—1); if there 
are the same number of particles, set rii — 0. Note that the 
radius tq was set to 8(7. Smaller as well as larger radii tq de- 
crease the performance of the order parameter as we found by 
comparing the order parameter with the actual order observed 
in the system. 
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Fig. 2 Sketch of the polar and equatorial regions used in the 
definition of the order parameter. 



Representative examples for the behavior of the resulting 
order parameter at constant angular frequency but increasing 
field strength are given in Fig. |3] As can be seen, in all 
the cases the value of ip grows with the field strength until 
it almost reaches a value of 1. Since the layers are usually 
not perfectly defined in our Brownian dynamics simulations, 
the order parameter typically takes on values that are slightly 
smaller than 1 even at very high field strengths. 

One also finds from Fig. [3] that there is a qualitative differ- 
ence in the behavior of t/j at high and low frequencies: The 
order parameter grows much more steeply at large frequen- 
cies, which means that the layers do not slowly emerge upon 
increasing the field strength but appear very rapidly. 



1.00 




Fig. 3 The order parameter ip at constant angular frequency (a) 
ujo = I, (b) 20, (c) 30, (d) 40. 

By inspecting snapshots corresponding to a given value of 
the order parameter, it turned out that the value i/jq w 0.6 may 
serve as an (approximate) lower limit for layer formation. 

Based on that criterion, we have scanned a broad range of 
frequencies and field strengths for the occurrence of layers. 
The results of this exploration of the parameter space are sum- 
marized in Fig. |4] Note that every simulation was started from 
a random configuration to avoid hysteresis-like effects. 

Within the layered "state", the translational structure of one 
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Fig. 4 Occurrence of layers depending on field strength and 
frequency of the driving field. The system parameters are chosen as 
described in Sec.|2] 



layer is disordered and becomes more and more homogeneous 
with larger ujq. In particular there is no pronounced hexagonal 
order as observed in earlier studies, ^ even though the particles 
tend to have six nearest neighbors at high ujq. This absence of 
pronounced in-plane order is probably a consequence of the 
low density considered (p* =0.1) and the Brownian random 
forces. Furthermore, depending on the initial conditions, we 
typically observe two or three layers in our simulation box 
(N = 864), which corresponds to an average vertical distance 
between the layers of about seven to ten particle diameters. 

The figure shows that the Uq — Bq diagram is separated into 
a layered and a non-layered region. Upon increasing the fre- 
quency from zero, the boundary first remains at roughly con- 
stant field strength, until it begins to rise with the frequency. 
This behavior is mirrored in Fig. |3] The larger the frequency, 
the higher the field strength at which the order parameter at- 
tains large values. 

A similar picture emerges from Fig. |5] where we have 
plotted the normalized absolute value of the magnetization, 
i.e. M/Mq = {\M{t)\)/Mo with M{t) = and 
Mq = Nfi, as a function of Bq for several frequencies. Note 
that |M(t)| is essentially independent of time for the states 
considered in Fig. |5] Clearly, the magnetization behaves dif- 
ferently in the regimes of small and large frequencies. One 
also finds from Fig. |5]that a degree of magnetization of more 
than 80 percent is required for layer formation to occur. 

In the following subsections, we will discuss the emergence 
and breakdown of layering in the different frequency regimes 
in more detail. Before doing so, it is worth to briefly comment 
on a technical issue encountered in our exploration of the pa- 
rameter space (see Fig. |4]i that concerns the behavior of the 



rotational temperature Trot = l/(2(-/V — 1)) J2^=i ^^i- Upon 
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Fig. 5 Absolute value of the magnetization normalized with respect 
to its saturation value over field strength at different rotational 
frequencies. The dots indicate after which point the system is 
considered layered according to our order parameter. 



More precisely, we define / as 
1 / 

-9(0,- (n+l)A(/.)y (14) 

where O is the Heaviside function, is the interval length 
to which we want to resolve the distribution, n is a positive 
integer or zero that satisfies nA0 < < (n + 1)A0, and, as 
before, (• • • ) denotes a time-average. 

We start by considering systems that are driven by fields of 
considerable strength (Bg = 10) with frequencies that admit 
layer formation (c/ Fig. |4]l. Results for the distribution / at 
three such frequencies luq are given in Fig. |6] For each value 
of cjq one observes a single, pronounced peak, reflecting a 
synchronized "state", in which the particles follow the field at 
constant phase difference. Note that the larger Wq, the larger 



increasing the diiving frequency loq from zero (at fixed B^), 
we typically also find Trot to increase, while its translational 
counterpart Ttrans ~ l/{i{N — 1)) Y^^=i '^^f stays approx- 
imately constant (close to the input value T). Similar tem- 
perature drifts have been observed in other non-equilibrium 
systems such as fluids in shear flow. In the latter context, 
the temperature is often redefined with respect to the differ- 
ences between the actual velocity of the particle and that of 
the flow field. Using a similar definition here (involving the 
difference between LOi and cjq), we find that this temperature 
is still not equal to T, but remains essentially constant over a 
broad range of frequencies. We also note that both the temper- 
ature drift and the actual location of the layer boundary in the 
Wq — Sq diagram depend on the chosen value of the rotational 
friction constant. 




Fig. 6 Distribution of the phase differences at Bq =10 and three 
frequencies uj'^ with = 7r/48. The systems are in layered states 
(see Fig.|4ll. 



3.3 Rotational dynamics in the layered state 

As mentioned earlier, the key argument for the appearance of 
the layers is that the time-averaged interaction between two 
fully synchronized rotating dipoles favors an in-plane configu- 
ration. In the following, we will investigate in more detail to 
what extent this assumption is actually fulfilled within the lay- 
ered region indicated in Fig. |4] To this end, we consider the 
distribution / of the phase differences 0; between the dipolar 
vector of particle i in the x — j/-plane and the external field. 



the phase difference between the particles and the field. This 
is not too surprising since an increase in the diiving frequency 
implies an increase in the rotational friction due to the (im- 
plicit) solvent and the presence of neighboiing paiticles. Fur- 
ther note that even though eqn (fl4l i contains a time-average, 
the phase distributions of these layered systems are essentially 
independent of time. 

To investigate the degree to which the particles actually ro- 
tate in the plane of the field, we also consider the distribution 
of the z-components of the angular frequencies 
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Fig. 7 Distributions of the z-component of the angular frequencies. 
Parameters are as in Fig. |6] The interval Auj is set to 1. 

In an ideal situation, in which the dipoles rotate perfectly 
with the field, the distribution g would have a single, sharp 
peak at a;* = ujq. Simulation results for g in the true many- 
particle system are shown in Fig. |7] where we have picked 
out the "states" already considered in Fig. |6] As expected in 
the layered regime, the functions g are characterized by one 
central peak around w cjq. However, we also see that there 
is a significant broadness in the distribution (as there is in the 
corresponding peaks of /). 

Finally, above a certain frequency, the layers disappear 
This is reflected in the emergence of a double-peaked struc- 
ture in the distribution of the phase differences, as illustrated 
in Fig. [8^. Moreover, we found that the non-averaged distri- 
bution of the phase differences is not independent of the time 
anymore. However, since we could not identify any system- 
atic time-dependence in this regime, we restrict ourselves to 
considering the averaged distribution. The first peak in / at 
(p « 7r/4 is due to particles that can still temporarily follow 
the field, whereas particles that are not able to do so anymore 
cause the structure of the rest of the distribution. The break- 
down of layering is also visible in the distribution g. Contrary 
to what is seen in a layered system, the angular frequencies of 
the majority of the particles are distributed around w as 
shown in Fig. [HJi. The much smaller peak at approximately the 
frequency of the external drive shows that only a small frac- 
tion of the particles follow the field at any given time. This 
fraction is further decreased as the frequency luq of the driv- 
ing field increases. Typical distributions / and g at values of 
ojq outside the layered regime are shown in Figs. and|9]5, 
respectively. Note that the roughly symmetric distribution of 
w* around approximately zero in Fig. [gji indicates that the 
particles are as likely to rotate in the direction of the field as 
they are to rotate in the opposite direction. 




Fig. 8 (a) Distribution of the phase differences of the system at 
Bq = 10 and uiq = 55.8, i.e. just outside of the region of layer 
formation. The resolution is — n/40. (b) Distribution of lj^ for 
a system at Bg = 10 and luq — 55.8 with Au — 1. 



Further note that at the large values of Bq considered in this 
section, the transition between states with the particles follow- 
ing the field at fixed phase difference and states where this is 
not true anymore happens in a very small range of frequencies . 

3.3.1 Effective single-particle theory. To understand the 
character of the high-frequency boundary between layered and 
non-layered states in more detail, we now aim to construct an 
effective theory that describes a single dipolar particle rotat- 
ing in a viscous medium. A similar consideration has been 
suggested for optically torqued nanorods by Shelton et al..— 
Clearly, such a single-paiticle approach cannot grant us direct 
insight into the formation of layers. However, it may help us to 
improve our understanding of the rotational dynamics isolated 
from many-particle effects. For simplicity, we assume that the 
rotational motion of the particle is restricted to the plane of 
the field and that it experiences rotational friction with fric- 
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Fig. 9 (a) Distribution of phase differences of the system at 

Bq — 10 and ujq = 60. The system is unlayered. (b) Distribution of 

oj* of the system at Bq = 10 and ujq = 60. The system is 

unlayered. 



tion constant 7. Then the rotational equation of motion is 

la — —7a + /ii?o sin(ajoi — a), (15) 

where a is the angle between the dipolar orientation and an 
arbitrary axis within the plane of the field. Equation (fTSt can 
be rewritten in terms of the phase difference i/i = wgi — a as 

/0 + 70 = 7^0 — /iSo sin(/). (16) 

We first consider the simplified case of negligible moments of 
inertia, i.e. an overdamped situation. Then eqn (fTST i reduces 
to the first order equation 

d(/) 
d7 



Wo 



= sill ( 



(17) 



where Wc = /ii3o/7 and r = ujci. This nonlinear differen- 
tial equation appears in various contexts such as the descrip- 
tion of overdamped pendula, superconducting Josephson junc- 
tions, and the synchronized emission of light by fireflies.^i^S 



For < Wo < "^c it has two fixed points characterized by 
(/) = (i.e. constant phase difference): One solution is a 
global attractor with <^ = arcsin(wo/a;c), and the other one 
is unstable with = tt — arcsin(ti;o/wc). These two solutions 
correspond to the phase differences at which the torque due to 
friction equals the torque that is due to the field. At wq = 
i.e. at (/) = 7r/2, the two solutions form a saddle-node bifur- 
cation and there are no fixed points for > Wc. At these 
high frequencies, the maximal torque that can be exerted by 
the field is insufficient to balance the frictional torque. The 
solution emerging after the bifurcation is a limit cycle with 
(i > 0. 
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Fig. 10 (a) The solid line shows the critical frequencies 
cjc = P^B^I^ji that are predicted by the single-particle theory in the 
BD frequency-field strength diagram (see Fig.|4). (b) Influence of 
the moment of inertia on the end of layer formation (dashed line: 
/* = 0.025, dotted line: 7* = 0.01). The solid line indicates the 
frequencies lj*. 

To which extent does the single-particle approach describe 
the true many-particle system of our BD simulations? In Fig. 
[TOk . the frequencies Wc (with 7 = see eqn (|5]l) are plotted 



This journal is ©The Royal Society of Chemistry [year] Jouma/ Nome, 2010, [vol], 1-FiTI |7 



into the Wq — Bp-state diagram (Fig. |4]i. At large frequencies 
loq and field strengths Bq, the straight line representing ojc has 
a slope similar to that of the boundary of the layered regime. 
This supports the idea that it is the (rotational) friction which 
eventually yields a breakdown of the synchronous rotations, 
and thus, the layering. A further observation from Fig. [TOi is 
that the true boundary frequencies (at given Bq) are somewhat 
smaller than uJc- One seemingly obvious reason for these de- 
viations is that the effective theory neglects any many-particle 
effects. Moreover, it does not take the Brownian random con- 
tributions into account that mimic the solvent "kicks" in eqn 
Both these factors could introduce perturbations of the ef- 
fective field that acts on a particle. Thereby the synchronized 
state could be destabilized already at frequencies uj < Uc- 
However, as it turns out, the more significant reason for the 
premature stop of layering is that the BD equations of motion 
involve (rotational) inertial terms, which are neglected in our 
single-particle approach. 

To check this point, we have performed additional BD sim- 
ulations with a lower moment of inertia (/* = 0.01). The 
resulting frequencies characterizing the boundary of the lay- 
ered state are shown in Fig. [Tob along with the original result 
(/* = 0.025) and the Une uJc- Clearly, decreasing the moment 
of inertia moves the true boundary substantially closer to the 
single-particle result. 

Finally, we note that the influence of the inertial (rotational) 
term can also be captured within our effective single particle 
theory. For / 7^ 0, eqn ( fTSI l can be written as 



dT'2 



dr' 



Wo 



(18) 



with v = j/y/JtBoI and r' = y/JiBoJlt. Similar to (fTSI l. this 
differential equation has a bifurcation at ujc, which means 
that the location of the line luc in Fig. [TOh remains un- 
changed. As before, the only stable solution at driving fre- 
quencies that are larger than ujc is a limit cycle. But addition- 
ally, for sufficiently small ly, it has a second bifurcation for 
some uj' with uj' < ujc as shown by Argentina et al. while 
investigating the transition between annihilation and preserva- 
tion of colliding waves.— This second bifurcation introduces 
a regime in which the limit cycle can coexist with the stable 
rotation. From the perspective of a many-particle system, one 
may speculate that the presence of the second solution per- 
turbs the rotation with constant phase difference (i.e., (p = Q). 

3.4 A density functional approach to layering in a per- 
fectly synchronized system 

We now consider systems at relatively low driving frequen- 
cies (ujq < 30), where, for sufficiently large field strengths 
Bq, the dipole vectors can follow the field in a perfectly syn- 
chronized fashion (see the discussion in the preceding sec- 



tion). According to our "phase" diagram in Fig. g] the field 
strength required to induce such synchronous and, at the same 
time, layered states, is about Bq « 4 — 6 for Wq < 30. The 
corresponding dipole-field coupling parameter ^BQ/k-oT = 
^*Bq/T* « 12 is significantly larger than the dipole-dipole 
coupling parameter (A « 6.7). Nevertheless, as seen in Figs. 
[3^ and b as well as Fig. |5l increasing Bq from zero at low 
driving frequencies yields a rather slow increase of the order 
parameter ip and the magnetization amplitude. 

Given the apparent interconnectedness between the rota- 
tional dynamics of the individual dipoles and the layering of 
the particles, we ask in the present section whether synchro- 
nization leads automatically to layering. Indeed, even in a 
perfectly rotating system, one would expect that the spatial 
symmetry-breaking associated with layering yields a decrease 
of translational entropy and thus may be unfavorable. 

To investigate this question we employ equilibrium density 
functional theory (DFT) for a system in which the dipole ro- 
tations are perfectly synchronized. Under such conditions the 
particles effectively interact via the time-averaged (inverted) 
dipolar potential given in eqn (fT2l i. By using this potential, the 
problem thus reduces to searching for an equilibrium phase 
transition in a system with effectively static interactions. 

Our density functional approach is based on the perturba- 
tion expansion of the free energy originally proposed by Ra- 
makrishnan and Yussouff in the context of fluid-solid transi- 
tions.— Up to second order in the density, the difference be- 
tween the Helmholtz free energy of a volume F of a system 
with non-uniform density p(r) and a reference system with 
homogeneous density po is given by— 



AT 
~V~ 



1 



dV[log(AV(r)) - 1] 
^ |^dVpo[log(A3po)-l] 



~ ^ /^'I'^i /'l'^2c(ri - r2)|p„ Ap(ri)Ap(r2). (19) 

In eqn (dUl, Ap(r) = p(r) - pQ with d^rAp(r) = 0, A 
is the thermal wavelength, and c(r)|po is the direct correlation 
function of the homogeneous system. 

Here we employ the random phase approximation (RPA) to 
calculate the direct correlation function. Assuming a hard 
sphere interaction in addition to the inverse dipolar potential 
JJ™ (eqn ^1), flie RPA amounts to setting 



c(r) = 



(r), r < a 
;3[/™(r), r>cr, 



(20) 



„PY 32 



where we used the Percus-Yevick direct correlation function, 
for the hard-sphere part. Note that within the RPA, the 
effects of the contributions of the long-ranged inverse dipo- 
lar interaction are treated in a mean-field fashion. To check 



8 I Journal Name, 2010, [voll, 14111 This journal is ©The Royal Society of Chemistry [year] 



this point, we have also calculated c(r) numerically by solv- 
ing the mean-spherical (MSA) integral equations. However, 
the changes in the free energies were found to be marginal. 

As a simple ansatz for the density profile in the layered 
state, we use 



p{r) = p{z) = po + pcos(fcz). 
Inserting this ansatz into eqn (fT9] l. we find 



(21) 



A 



1 



^\zp{z)\og(P^ 
\ Pa 



>^L^m. (22) 



where /S.F is the free energy of the volume A • Al, ^ is an 
area in the x — y-direction and — 27r/fc. Further, c is the 
Fourier transform of c and c{k) = c{kez). In the RPA, we 
have 



£(fc)=47r( / drr2jo(fcr)cPY(r)+/i2/3:^i^ 



, (23) 



where j„ are spherical Bessel functions of order n. (For the 
treatment of the dipolar interactions in eqn (l2Jt . see Ref. [33h 
We now use eqn (|22] | to search for a phase transition be- 
tween the homogeneous and the layered state. In principle, 
this search requires a minimization of AJ^/A with respect 
to both the parameters p and k that characterize the inhomo- 
geneity of the system (see eqn (|2T])). It turns out, however, 
that AF /A becomes minimal with respect to k for fc — >■ 0, 
which corresponds to an infinite distance between the layers. 
Clearly, this is not compatible with the implicit assumption 
of a finite wavelength. Therefore we have fixed the parame- 
ter k = 271 /\l to physically reasonable values, i.e. to val- 
ues suggested by our BD simulations. At p^ = 0.1, we find 
an average layer distance of approximately 7.2cr (see below). 
This leaves the coefficient p as the only minimization param- 
eter. Results for the function AJ^(p)/(j4Ai) with fixed dis- 
tance \l = 7.2(7 between the layers at various values of the 
parameter /i* are plotted in Fig. [TTh . 

The different curves in Fig. [TTk reveal a behavior typical 
of a second-order phase transition. For p* < 2.27, the free 
energy has only one minimum at p = corresponding to an 
homogeneous state. This changes at p* « 2.27: For larger 
values of p*, the solution at p = represents a maximum, 
and the only minimum occurs for p > 0. The corresponding 
negative values of AT /A indicate that it is indeed the layered 
state which is now globally stable. 

We have repeated the DFT calculations for a number of den- 
sities in the range 0.01 < Pq < 0.4. To find reasonable values 
for the corresponding wavelengths Al in the layered state, we 
ran BD simulations at fixed dipole moment p* — 3.4, fre- 
quency luq = 8, and field strength Bq — 50. With this choice 
of the parameters, the particles are almost perfectly aligned, 
justifying the key assumption of our DFT approach. Fitting 
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Fig. 11 (a) Free energy difference between a layered and a 
homogeneous state as a function of the coefficient p (see eqn \2li ) 
at different values of the parameter p* . The case p = corresponds 
to the homogeneous solution. The overall density is set to po = 0.1. 
(b) Equilibrium phase diagram of a perfectly synchronized system. 
The gray area indicates the stability range of the layered state 
according to our DFT calculations. Also shown are BD results (at 
cjq = S, Bq — 50) with the open circles (solid triangles) 
representing layered (non-layered) states. 



the resulting distances as functions of po, we found the ap- 
proximate relation d/cr 1.05p*~"'^*, which was then used 
as an input in the DFT (i.e., Al = d). 

The resulting phase diagram in the pp — p* -plane is plot- 
ted in Fig. [TTb . It is seen that the DFT predicts a layering 
transition for all but the smallest densities (pq > 0.01) in the 
shown parameter range, with the actual values of p* varying 
substantially with p^. Indeed, the lowest threshold is found 
at Pq « 0.2. Also shown in Fig. [TTb are BD results for the 
appearance of layers in nearly perfectly synchronized systems 



(wq — 8, — 50) at various values of p*. As in Sec. 13. 21 the 
presence of layers was detected on the basis of the order pa- 
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rameter defined in eqn (fTsl i. yet with a slightly different defini- 
tion of the cutoff radius entering the order parameter (tq = d). 
Comparing BD and DFT, it is seen that the DFT predicts the 
true phase boundary in perfectly synchronized systems in a 
qualitatively correct manner (including the strong increase of 
/X* upon Pq — > 0). Moreover, the DFT results are also quite 
reasonable from a quantitative point of view. 

From a physical perspective, clearly the most important 
conclusion is that even in a perfectly synchronized system, 
a sufficient decrease of interaction energy stemming from the 
time-averaged dipolar potential is required to overcome the 
entropy cost due to layering. 

Finally, we briefly discuss our DFT results in the light of 
a recent Monte Carlo study by Smallenburg and Dijkstra, ^' 
who obtained full equilibrium phase diagrams of systems in- 
teracting with inverted dipolar interactions. To model the 
short-range part of the interaction, Smallenburg and Dijkstra 
used either just hard spheres or hard spheres with an addi- 
tional Yukawa repulsion.^' In the first case, layer-like struc- 
tures were only observed in the gas-liquid coexistence region. 
On the contrary, the Yukawa system exhibits a stable layered 
phase with fluid-like in-plane structure. Comparing these lat- 
ter results to our DFT predictions, we find that the predicted 
strength of the inverted dipolar interactions required for layer 
formation is indeed comparable. On the other hand, we find 
the onset of layer formation at much lower densities. Apart 
from the obvious approximations in our theory, we also at- 
tribute these deviations to the fact that the repulsive Yukawa 
interaction used in Ref. QAj is much stronger than our soft 
sphere one. 

4 Conclusions 

In this study we have combined BD computer simulations, an 
effective single-particle theory, and an (equilibrium) density 
functional approach to explore the dynamic behavior of sys- 
tems of dipolar particles in planar rotating fields. 

One main result from our BD simulations is a non- 
equilibrium "phase" diagram identifying the domain of lay- 
ered states in the ujq-Bq plane (at constant particle density 
and dipolar coupling strength). At low driving frequencies, 
the change from unlayered to layered (and fully synchro- 
nized) structures occurring upon increase of Bq is related to 
a quasi-equilibrium phase transition, i.e. a many-particle phe- 
nomenon. The transition is induced by the competition be- 
tween the time-averaged, inverted dipolar interactions favor- 
ing in-plane configurations and the loss of translational en- 
tropy accompanying the one-dimensional translational order 
While this competition also occurs for systems of polarizable 
or superparamagnetic particles, the additional complication in 
the present system of permanent dipoles is that the field first 
needs to overcome the dipolar fluctuations. Though we have 



neglected this issue in our DFT approach, we would expect 
that the fluctuations just shift the transition predicted by the 
DFT towards larger field strength. 

Completely different behavior is found at high frequencies 
and field strengths. Under these conditions, the picture of syn- 
chronously rotating dipoles (with constant phase difference 
relative to the field) breaks down. Instead, one observes a mix- 
ture of rotating and counter-rotating or resting particles , as our 
analysis of various angle distributions reveals. The desynchro- 
nization induces, at the same time, a breakdown of the transla- 
tional, layered structure. Despite this complex many-particle 
behavior, we have shown that the boundary can be well de- 
scribed in terms of the critical frequency Wc(i?o) that arises 
from a bifurcation in an effective single-particle approach for 
the rotational motion in a viscous medium. This indicates that 
the appearance of the high-frequency boundary is essentially 
a friction-induced effect. 

A similar frequency-induced desynchronization effect has 
recently been discussed by Hartel et al.,— who investigated a 
system of interacting elongated particles in a rotating electric 
field via dynamic density functional theory. Assuming a con- 
stant number density, the important dynamic variable within 
the density functional approach is the orientational distribu- 
tion as function of time. At low and very high frequencies, the 
distribution behaves similar to our distribution / in that there 
is either a single peak (reflecting synchronized motion with a 
constant phase difference) or no peak at all. In the transition 
regime, however, Hartel et al. detected various new dynamic 
states characterized by time-dependent oscillations and split- 
ting of the peak in the distribution as well as an overtaking 
by the driving field. In the present study we did not observe 
such states, not even when looking at the time-dependence of 
our orientational distributions (or the magnetization). It re- 
mains to be investigated whether these qualitative differences 
in the rotational motion of anisotropic many-particles systems 
are just due to differences in the specific model system, or 
due to the fact that our results are based on a microscopic ap- 
proach rather than on the density field approach used in Ref. 
|23[ Indeed, the relation between the microscopic and meso- 
scopic dynamics in driven systems is an issue also discussed 
in other, related contexts, such as the shear-induced dynamics 
of nanorods. ^'^ 

Finally, it is worth to briefly comment on the relevance of 
our dimensionless model parameters in the context of real sys- 
tems. The equilibrium parameters considered here (density 
p* = 0.1, dipolar coupling strength A w 6.7) correspond to 
those of a strongly coupled ferrofluid exhibiting chain for- 
mation. Regarding the driving field, however, most of our 
dimensionless frequencies Wq are probably beyond the cur- 
rently accessible range. In many experiments involving rotat- 
ing fields, the size of the (typically superparamagnetic) par- 
ticles considered is about 1 /im.^^ A driving frequency of 
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Wq = 10 (which is well within the layered domain) then cor- 
responds to an actual frequency of about 10 kHz if we assume 
room temperature (T = 293 K) and a mass density of 5 g/cm^. 
This is 1-2 orders of magnitudes larger than the frequencies 
used in the literature/ " FerrocoUoidal particles, which have 
permanent dipoles (such as the ones considered here), are of- 
ten much smaller with sizes of about 10 nm. In that case, 
Wq = 10 corresponds to a driving frequency of about 1 GHz. 

These considerations suggest that realistic driven systems 
will be fully synchronized and layered according to our 
"phase" diagram in Fig. |4] We note, however, that the ac- 
tual location of the desynchronization line encountered upon 
increasing luq depends on the friction constant used in our BD 
simulations; i.e., increasing the friction constant shifts the line 
towards lower frequencies (consistent with the single-particle 
theory). Moreover, we have neglected in our study the many- 
particle character of the hydrodynamic interactions induced 
by the solvent. We would expect these interactions to effec- 
tively increase the friction and thus shift the boundary towards 
even lower frequencies. Clearly, it would be very interesting 
to actually incorporate such interactions by using refined sim- 
ulation methods such as, e.g., stochastic rotation dynamics.-'^ 
Hydrodynamic interactions may also be relevant to better ex- 
plore phenomena such as chain-to-cluster transitions that have 
been revealed by recent studies."* These issues, as well as the 
dynamic behavior in even more complex field geometries, will 
be the subject of future studies. 
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